Equation of State of nuclear matter in a Virial expansion of nucleons and nuclei 
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We study the equation of state (EOS) of nuclear matter at subnuclear density in a Virial expansion 
for a nonideal gas. The gas consists of neutrons, protons, alpha particles, and 8980 species of nuclei 
with A > 12 and masses from the finite-range droplet model (FRDM). At very low density, the 
Virial expansion reduces to nuclear statistical equilibrium. At higher density, the Virial results 
match smoothly to the relativistic mean field results discussed in our previous paper. We tabulate 
the resulting EOS at over 73,000 grid points in the temperature range T — 0.158 to 15.8 MeV, the 
density range ub = 10~* to 0.1 fm~^, and the proton fraction range Yp — 0.05 to 0.56. In the 
future we plan to match these low density results to our earlier high density mean field results, and 
generate a full EOS table for use in supernova and neutron star merger simulations. This Virial 
EOS is exact in the low density limit. 
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I. INTRODUCTION 

One of the main ingredients in simulations [H of 
core coUpase supernovae and neutron star mergers is the 
Equation of State (EOS) for hot dense nuclear matter. 
The EOS, along with detailed information on the com- 
position of nuclear matter, may play an important role 
in the neutrino-matter dynamics [3j and supernova ex- 
plosion mechanism. In a previous paper |j| we used 
a relativistic mean field (RMF) model to calculate the 
free-energy of non-uniform matter at intermediate den- 
sity and uniform matter at high density. 

Simulations of the core collapse and explosion of su- 
pernovae depend heavily on the EOS at subnuclear den- 
sity, especially the detailed composition. In this work, 
we study subnuclear density nuclear matter in a Virial 
expansion for a nonideal gas, consisting of neutrons, pro- 
tons, alpha particles, and 8980 species of heavy nuclei 
{A > 12) with masses from the finite-range droplet model 
(FRDM) 0. The Virial results will cover the density 
range ub = 10~^ to 0.1 fm~'^, the temperature range T 
= 0.158 to 15.8 MeV, and the proton fraction range Yp 
= 0.05 to 0.56. (For temperature higher than 15.8 MeV, 
matter is uniform and fully described in the RMF model 
[3|-) The distribution of nuclei for given conditions is ob- 
tained in this approach, while the existing EOS tables of 
Lattimer-Swesty (L-S)j6iand H. Shen, Toki, Oyamatsu 
and Sumiyoshi (S-S) 0,111, ^ single heavy nucleus 
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approximation. Our Virial EOS will be matched, using 
a thermodynamically consistent interpolation scheme, to 
the RMF EOS obtained in our previous paper Q to gen- 
erate a full EOS table for supernova. 

Detailed information on the distributions of nuclei in 
the EOS table is important for neutrino-matter dynam- 
ics. Neutrinos radiate 99% of the energy released in su- 
pernovae. Besides gravitational wave signals, neutrinos 
are the only messenger through which one can directly 
probe the EOS inside supernovae. The neutrinosphere is 
the surface of last scattering before ly or P escape. The 
neutrinosphere is expected to occur at a density around 
10^^ g/cm'^ and this is consistent with the available infor- 
mation from a handful events in SN1987a The com- 
position of matter at subnuclear density constrains the 
position of the neutrinosphere and infiuences the spectra 
of emitted neutrinos and antineutrinos. For example, in 
a recent study [l3|, fight nuclei with mass 2, 3 and 4 were 
found to have an important influence on the spectra of 
anti-electron neutrinos. 

For matter at subnuclear density and low entropy, 
heavy nuclei tend to form. Nuclear statistical equilib- 
rium (NSE) models treat low density nuclear matter as a 
system of noninteracting nuclei in statistical equilibrium, 
taking into account the binding properties of heavy nu- 
clei. This has been widely used in nuclear astrophysics 
. Recently, there have been several NSE based studies 
of the supernova EOS, see for example Refs. [1^ and 
These studies use modern mass tables with thousands 
of nuclei and include excluded volume effects [13]. The 
NSE models have the advantage of generating thermody- 
namically consistent EOS tables. However they also have 
several disadvantages. First, NSE models themselves can 
not be used to describe well the non-uniform matter at 
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nearly nuclear density, when exotic pasta phases may ap- 
pear. To generate a complete EOS table, NSE models 
need to be matched to uniform matter models at high 
density. Thus low and high density matter are usually 
described with different models. Second, NSE models do 
not provide a systematic way to include the strong in- 
teractions between nuclei. These interactions could be 
important as the density increases. Moreover, neutron 
matter at low density is closely analogous to a unitary 
gas, whose properties can not be described satisfactorily 
in either NSE models or normal mean field approaches 
(see for example The Virial expansion, on the other 

hand, has been successfully used to describe the proper- 
ties of unitary gases [11] as well as neutron matter at low 
density [lB|- 

In a previous work [l7| on low density nuclear matter 
consisting of neutrons, protons and alpha (a) particles, 
the Virial expansion has been used to systematically in- 
corporate second virial corrections from nucleon-nucleon 
(N-N), Na and act elastic scattering. These light nu- 
clei components are expected to dominate in matter with 
high entropy. The effects of second virial corrections are 
found to be important for the composition and for the 
equation of state. In this work we use the Virial expan- 
sion for a nonideal gas consisting of neutrons, protons, al- 
pha particles, and thousands of heavy nuclei. The heavy 
nuclei included in this work are from FRDM mass tables 
0, which have 8979 nuclei with A > 16. We also include 
"'^^C in the mass table. As mentioned above, the mass 2 
and 3 nuclei, and other A < Yl nuclei asside from ^He 
may be important for the spectra of anti-electron neu- 
trinos. The inclusion of these light elements will require 
information on second virial corrections among them and 
nucleons, which are under current investigation [l^. In 
this work the light elements are represented by alpha par- 
ticles, which typically have the largest abundance among 
the light elements in statistical equilibrium. At very low 
density, the Virial expansion agrees with nuclear statisti- 
cal equilibrium models. As the density rises, we find that 
the Virial gas matches smoothly to previous mean field 
results. 

In this work we include the second order virial correc- 
tions among nucleons and alphas as in Ref. 17]. Par- 
tition functions for heavy nuclei are included using the 
recipes of Fowler tt aZ.[l^. In addition, some calcula- 
tions are presented using partition functions based on 
the recipe of Rauscher et aZ.fl^ for comparison. Equiv- 
alently, the partition functions can be considered as the 
sum of successive high orders of the Virial expansion for 
heavy nuclei. There have been many studies on the level 
density and partition functions of hot nuclei in astrophys- 
ical environments. For large scale astrophysical applica- 
tions, it is necessary to find both reliable and computa- 
tionally practical methods for the level density. Most of 
these studies p^|-[23l] followed the orig inal non-interacting 
Fermi gas model of Hans Bethe [2J]. For astrophysical 
nuclear reactions with temperature below a few times 
10^°X [13, d^, this phenomenological approach gave ex- 



cellent agreement with more sophisticated Monte Carlo 
shell model calculations [2^ , as well as combinatorial ap- 
proaches [H, [23] ■ This justifies the application of the 
Fermi-gas description at and above the neutron separa- 
tion energy. For temperatures higher than a few times 
lO^^X, there are big ambiguities in the values of the par- 
tition functions. However, as suggested by some authors 
12] and supported by our own calculations, these uncer- 
tainties have only a small impact on the thermodynamics 
of dense matter. 

The effects of Coulomb interactions in the plasma can 
be estimated by the plasma parameter Fp = [Ze)"^ /akT, 
where Z is the atomic number of the nucleus, T is the 
temperature and a is the spacing between nuclei. For 
matter at low density, Tp is smaller than one and the ef- 
fect of Coulomb corrections is small. However for matter 
at higher density and when the dominant species carry 
large charges, Tp can be much greater than one and the 
effect of Coulomb interactions should be taken into ac- 
count. The Coulomb correction to the plasma has been 
studied analytically up to high Tp by the cluster expan- 
sion [2^ . Generally the correction due to electron- ion in- 
teractions will reduce the free energy of the plasma and 
eventually crystalize the matter at high density. For sim- 
plicity, in this work the Coulomb interactions between 
nuclei and electrons are included via a Wigner-Seitz ap- 
proximation with effective ion spheres for each species of 
nuclei, wherein local electrical neutrality is maintained. 
This Wigner-Seitz approximation for the Coulomb cor- 
rection will be compared with a more rigorous cluster 
expansion method. 

Based on the above Virial expansion, we generate an 
equation of state table which covers the range of temper- 
atures, densities and proton fractions shown in Tab. ID 



TABLE I: Range of temperatures, densities and proton frac- 
tions in tlie EoS table. 



Parameter 


minimum 


maximum 


number of grids 


T [MeV] 




10^'^ 


20 


logio(ns) [fm^^] 


-8.0 


-1.0 


71 


Yp 


0.05 


0.56 


52 



There are 7 points in temperature (0.16, 0.26, 0.40, 
0.63, 0.71, 0.79 and 0.89 MeV) for temperature below 1 
MeV. For higher temperatures we use a spacing of 0.1 
in logio(r/ [MeV]). This is a total of 20 different tem- 
peratures from T = 0.16 to 15.8 MeV. We use a spacing 
of 0.1 in logio("-B/ [fm '^]), giving a total of 71 points 
in density for ub = 10^^ 0.1 fm"'^. We use a linear 
step 0.01 in proton fraction, giving a total of 52 points in 
proton fraction for Yp = 0.05 to 0.56. There is a total 
of 73,840 data points in the Virial gas calculation. This 
took 1,000 CPU days in Indiana University's supercom- 
puter clusters. 

The paper is organized as follows: in section |ll] our 
Virial expansion for a nonideal gas is explained in de- 
tail. In section IIIII we present the recipes for the nuclear 
partition function used in the Virial expansion. Section 
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IIVI discusses the effect of different mass tables on the 
EOS, and shows several examples for the free energy and 
the distribution of nuclei in the Virial EOS. Section |V] 
presents a summary of our results and gives an outlook 
for future work. 



where ta = I.IGA^^^ fm is the nuclear radius (in accor- 
dance with that in FRDM ^]), and is the average ion 
sphere radius defined through 



(6) 



II. FORMALISM 

We now describe our Virial expansion formalism for a 
gas consisting of neutrons (n), protons (p), alpha parti- 
cles and heavy nuclei. The grand partition function Q for 
a gas at pressure P and volume V is expanded to second 
order in the neutron fugacity z„, the proton fugacity Zp 
and the alpha particle fugacity follows, 



logg 
V 



P 2 

7^ = Tjt^" + + {'4, + ^n)^n + 2ZpZj)pn] 



(1) 



where fli is the partition function for nuclei and 
bn , bpn , ba , and bnn are the second Virial coefficients as 
defined in Ref. hj. The sum on i runs over different 
heavy nuclei, for which we use the FRDM mass table Q 
for A> 16 and include ^^C. The thermal wavelength for 
species a is Aa, 



Aa = y/2T: /niaT, a 



p, a, nuclei. 



(2) 



From now on i, j, ... are used for sums over heavy nuclei 
and a, &, ... are used for sums over all species. 

There exist several recipes for the nuclear partition 
function fi^. We will use that of Fowler et a/. [13] in 
this work. We also consider the choice of Rauscher et 
aL[20| as an option. Different choices of partition func- 
tions change the matching densities to our mean field 
results slightly, but the infiuence on the thermodynamics 
is negligible. This is also discussed in Ref. [l^ . 

Chemical equilibrium between nucleons and a heavy 
nucleus with Z protons and N neutrons insures. 



(3) 



where fii, fip, fin are chemical potentials of the heavy nu- 
cleus, protons and neutrons, respectively. Therefore the 
fugacity of a heavy nucleus is readily obtained. 



exp(^i + Ei)/T 



(4) 



where Ei is the binding energy of heavy nucleus "^Z. 

We consider the Coulomb interaction between the elec- 
tron background and nucleus following Baym et aZ.[29j. 
but generalized to multiple species of nuclei. The total 
Coulomb energy of a nucleus in an electron background 
is, 



3 Zfa , 

5 TA 



3 

2 n 



(5) 



We emphasize that the sum over j runs over heavy nuclei. 
This definition ensures local charge neutrality inside each 
ion sphere. The first term in Eq. ^ comes from the 
proton-proton Coulomb energy in the nucleus and is also 
included in the binding energy of the nucleus. When 
ta approaches the total Coulomb energy in a heavy 
nucleus vanishes as expected. So the Coulomb correction 
to the binding energy of nucleus is 



5 Tyi 2 Ti 2 ri 



(7) 



Adding this correction to the binding energy of a nucleus, 
Eq. (H]) becomes 

= exp(M. + K,-£;f)/r = z^^z^et^--^?)/^. (8) 
The densities of each species can be obtained from 
/ d logQ 



\dZa V 



(9) 



This gives, 



2 

— "\3" '^^n^n ^~ '^^p^n^pn ~l" ^^a^n^an\f (10) 

1 



A3 



{Zn + Zp)b 

O 



(12) 
(13) 



The mass fraction of species a{A, Z) is defined as, 

Xa = AaUa/uB. (14) 

The two conditions used to determine the fugacities of the 
neutrons and protons are that the total baryon density 
is conserved, 

ns = n„ -I- np -f 4na + ^ AjUj, (15) 

i 

and that the proton fraction is given by, 

Yp = {np + 2na + ^^Zini)/nB. (16) 

i 

Since the Coulomb correction is included, one extra loop 
is needed to self-consistently determine the values of the 
ion sphere radii for each species. 
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The entropy density s of the Virial gas is obtained 
from, 



III. NUMERICAL DETAILS 



dP 
df 



5 P 2T 



pni 



T 



^2 IT'i^OgZi + n„logZ„ + UplogZp + UalogZa 

dE? . 



dT 



(17) 



where the prime indicates a derivative with respect to 
temperature. Note that the last term is the Coulomb 
correction to the entropy, which is nontrivial to evaluate 
directly. However the free energy of the Virial gas can 
be obtained directly (see discussion below) so that we 
can evaluate the entropy from the derivative of the free 
energy. 

The energy density e can be calculated from the en- 
tropy density, 



(18) 



where the sum on a runs over neutrons, protons, alphas 
and heavy nuclei. 

The free energy density / is given by 

/ = e - Ts = ^ na/ia - P 

a 

= n„r log Zn + ripT log Zp + n^T log Za - riaEa 
+ [n^T\ogz, - n,{E, - Ef')] - P. (19) 

i 

From the above equation we define an effective Coulomb 
correction to free energy per nucleon A//v4, 



riB 



The free energy per nucleon is 

F/A = f/uB. 



(20) 



(21) 



The thermodynamic pressure P can be obtained from 
the free energy, 



Pf 



th 



which can be rewritten as 



(22) 



(23) 



The second term in Eq. ((23)) is the correction to thermo- 
dynamic pressure due to the Coulomb interaction. We 
present results for the Virial equation of state in Section 

EYl 



Pth = P + ^'^BY^ni-^\T,Yp 



In Subsection llll Al we describe some details of the eval- 
uation of the partition functions and then in Subsection 
HUB I we briefly describe the parallel computations of the 
EOS for many different temperatures, densities, and pro- 
ton fractions. 



A. Recipes for partition function 

Fowler et a/.[l^ proposed an efficient approximation 
for the partition function of hot nuclei, which has a closed 
form. 



dEp{E)cxp{-E/T) - nc 



(24) 



where fid is the contribution from known discrete states 
and flc is the continuum subtraction. One could use 
fid — (2Jo + l) where Jq is the ground state spin. Inaccu- 
racies from this approximation will become progressively 
unimportant beyond 10^*^ K. The continuum subtraction, 
on the other hand, only becomes important beyond 10^^ 
K (see also [30]), by this temperature uniform matter will 
be more stable in most regions of phase space, as will be 
shown in the following discussions. So in the tempera- 
ture range we are interested, one can discard the latter 
term. Therefore the nuclear partition function becomes, 



(2Jo + l)+ / dEpiE)exp{~E/T), (25) 



E. 



and its derivative versus temperature is, 

cEt 



TO^ = / ' dEp{E)exp{^E/T)^. 

J Eh ^ 



(26) 



A widely used expression for the level density p{E) is the 
back-shifted Fermi gas formula [2^ . 



p{E) 



/tt cxp(2\/ aU) 
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l/4[/5/4 



(27) 



where a is the level spacing, and U = E — S with 6 a back- 
shift parameter related to pairing. Various prescriptions 
of these two parameters are available, which reproduce 
more rigorous results for the level density from Monte 
Carlo or combinatorial calculations. We will use the pre- 
scription from Fowler et al |19l| in our calculation and 
include that of Rauscher et aL[20| as an option in our 
code. 

The integral limits are determined by the following 
equations ^19.] : 



Ed - imin(5„,5,) 



(28) 



Et = mm{Sn + ER,Sp + ER + -E,). (29) 
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Sn and Sp are the single neutron and single proton sep- 
aration energies, whicli could be obtained from the mass 
table. Ell is the zero-point energy, E^i — /2ME? 
with R = l.2b{A - 1)1/3. The Coulomb barrier is 
Ec = {Z — \)a/R. For unstable nuclei {Sn or Sp < 0), 
the partition function is set to the ground state value fid- 
There is still a large ambiguity regarding the value of the 
integral upper limit Et-, which is related to the contribu- 
tion of the continuum. This ambiguity will introduce big 
changes in the partition function only when the temper- 
ature is above several MeV (few times lO^^if). We com- 
pared results with different choices of Et and found the 
ambiguity changes the transition density to our relativis- 
tic mean field EOS ^ slightly, but the influence on the 
thermodynamics is negligible. This was also discussed in 
a previous work [l"2j . 

When S is larger than Ed, 5 — Ed should be used. In 
this case Eq. (P7|) is inapplicable due to the pole in the 
lower integration limit. A constant temperature formula, 

piE) = C7exp(C//r,), (30) 

is used instead. The constants C and Tc are obtained by 
the continuity of the level density and its first derivative 
when matched to Eq. (P7| . 

1. Choices of a and 5 by Fowler et al.. fT^ 

In our calculation we use the prescription of Fowler et 
al.^W] for the level spacing parameter a (in MeV"^) and 
pairing parameter 6 (in MeV): 

Z<30: a ^ 0.052^1-2, S = dp-80/A, , . 
Z>30: a = 0.125A, S ^ Sp - 80/ A ~ 0.5, ^ ^ 

where Sp = (11 • A-i/2MeV) [l + i(-l)^ + i(-l)^] . 
When S > Ed, S is set to Ed, and Eq. ([50)1 is used below 
an energy 2Uc {Et > 2Uc > 0). By continuity of the level 
density and its derivative, the constants C and are 
obtained, 

1 _ 5 1 Vfl 

C = ^a'^f^U-^/^eMl + ^/^c). (32) 
In our calculation Uc is chosen to be S. 

2. Choices of a and 5 by Rauscher et al. . 

In the parameterization of Rauscher et al., the level 
spacing a depends on energy, 

1 - e-'''^ 

a{U,Z,N) = {aA + pA^/^)[l + b ], (33) 

where b = E^ic describes properties of a nucleus differing 
from the spherical macroscopic energy, which is already 



given in the FRDM mass table. Based on the FRDM 
mass formula, the best fitted values for a, /3, and 7 are 
obtained by comparing the level densities with experi- 
mental analysis: a = 0.1337, P = -0.06571, 7 = 0.04884 

The pairing parameter S is given by, 

S^^{An{Z,N) + Ap{Z,N)}, (34) 

where the neutron pairing energy is 

An{Z,N) = ^[2E^{Z,N)-E'^{Z,N-l)-E^{Z,N+l)] 

(35) 

and E^{Z,N) is the binding energy of nucleus {Z,N). 
The proton pairing energy Ap{Z,N) is calculated in a 
similar way. The constants C and Tc are obtained as 
following, 

Tc AUc \l Uc -i a y a ' 

C = p(C/,)e-^=/^% (36) 

where 

a' = (a^ + /3A2/3)A(_i + e-T^=+7i7,e-^^=). (37) 



B. Computational methodology 

We calculate, in parallel, the points in parameter space 
covered by the Virial EOS, analogous to the way used 
in our previous paper to perform relativistic mean field 
calculations 4] . There are a total of 73,840 points in the 
3-dimensional parameter space of temperature T, proton 
fraction Yp, and density ub- Each point takes up to 20 
minutes to calculate. The overall Virial EOS took about 
1,000 CPU days in Indiana University's supercomputer 
cluster. 

Each point in the parameter space was mapped to 
a unique integer that we refer to as the job index. A 
file, runlist, was prepared with a list of job indicies 
for the whole parameter space, and a single character 
(A=available, R=running, r=Re-running, C=complete, 
T=time-limited and F=failed) that gives the status of 
calculations for that job index. A Message Passing Inter- 
face (MPI) parallel wrapper code manages the running 
of the many requested tasks. Typically, one parallel job 
requests a set of compute cores (usually 256). Each MPI 
rank, using a single CPU core, is assigned one job index 
corresponding to one point in the parameter space and 
evaluates the required quantities. 

Initially, rank zero of the MPI job 

• locks the job listing file runlist, 

• reads runlist until a list of available tasks is filled, 

• closes runlist and releases the lock. 
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• passes a job index to each MPI rank and begins the 
calculation for that job index. 

When the calculation completes (or time-limits or fails) 
for a given MPI rank, the status character for the job in- 
dex in runlist is modified appropriately. The now avail- 
able MPI rank will search runlist for the next available 
task and the calculation restarts for the new job index. 
Since completion occurs asynchronously file locking is not 
used for this part of the process. 

A simple batch job runs through the points in param- 
eter space. A wall clock limit (48 hours) is used. Each 
rank of the MPI job can run a series of points using the 
above procedure, efficiently using each available core for 
the requested wall clock period. One job per core is run- 
ning when the wall clock limit is reached. These jobs are 
identified by being left in the "R" state after the batch 
job completes. This procedure allows us to calculate > 
99 % of the points in the runlist file. 



IV. RESULTS 

In this section we discuss the Virial EOS in detail. 
First we show the effect of different mass tables on the 
EOS and its influence on the matching to RMF results. 
We also discuss the coulomb correction and compare our 
results with some analytical cluster expansion analysis. 
Second we discuss the matching between the Virial EOS 
and RMF EOS for several choices of temperature and 
proton fraction. Finally we show some examples of the 
distribution of nuclei obtained in the Virial EOS. We 
also compare some examples of mass fractions of nuclear 



matter between Virial EOS and existing EOS tables. 



A. Effect of different mass tables 

We use both the FRDM Q and HFB14 [sH mass tables 
in the Virial expansion. Here we compare the free energy 
and average charge number Z (average mass number A = 
Z /Yp) from the Virial expansion with the two different 
mass tables. We also match them to the relativistic mean 
field results from our previous paper 



1. Free energy 

In Figure [TJ free energies are shown as a function 
of density for nuclear matter at T = 1 and 3.16 MeV 
and Yp = 0.3. Virial gas results are obtained from 
Eq. ^ for the FRDM (black solid curve) and HFB14 
(red dashed curve) mass tables. Nonuniform (or Hartree) 
mean field calculations (circles) and uniform matter cal- 
culations (dotted dashed curve) are also shown (See Eqs. 
(10) and (18) in Ref. |:4i]). These calculations give lower 
free energies at high densities. The two different mass 
tables in the Virial expansion give very similar free ener- 
gies for nuclear matter. The difference between the black 
curve (FRDM) and red curve (HFB14) is very small until 
densities where the Hartree mean field results are more 
stable than either mass formula. For these two cases, the 
transition from Virial to Hartree EOS occurs at a density 
of about 3.98x10^3 fj^-s^ 



J 
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FIG. 1: (color on line) Free energy per nucleon of nuclear matter at T = 1 (a) and 3.16 MeV (b) with Yp — 0.3. 

I 



2. Average A and Z 



EOS with either FRDM or HFB14 mass tables, and from 



Similar to Fig. [TJ the lower panels in Fig. [5] give the 
corresponding average charge number Z from the Virial 
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Hartree mean field results. Again the Virial EOS with 
the two mass tables predict very similar values for Z. 
Moreover, the Virial EOS with either mass table gives 
very similar Z to that from the Hartree mean field results 
at the transition density 3.98x10^"^ fm^'^ (blue dotted 
line). The fiuctuation of Z in Hartree results below the 
transition density (at T = 3.16 MeV) is probably due 
to finite step error in search of cell size. In conclusion, 
the composition, free energy, and transition density to 
Hartree mean field results, depend little on the choice of 
mass tables. 



In the two upper panels in Fig. [2l we compare 
the Coulomb energy correction in our Virial expansion, 
Eq. (I20p . with an analytic cluster expansion for the one- 
component plasma, Eq. (22) in '^s']. The overall agree- 
ment for densities below the Virial-Hartree transition is 
good, though at higher temperature the differences be- 
come larger. The reason is probably that we calculate the 
multi-component contribution to the Coulomb correction 
in the Virial gas while we only used the average charge 
number in the analytic formula for the one-component 
plasma. 




FIG. 2: (color on line) Upper panels (a), (c) show Coulomb corrections; lower pannels (b), (d), show average charge number of 
heavy nuclei, in nuclear matter at T = 1 MeV, left panels (a), (b), and r=3.16 MeV, right panels (c), (d). The proton fraction 
is Yp = 0.3. 



B. Free energy and phase boundaries 

We show in Fig. [3] the transition densities between the 
Virial and RMF EOSs. We also show the free energy 
per nucleon F jA as a function of density for T = 
1, 3.16, 6.31 and 10 MeV. At low densities, FjA is ob- 
tained from Eq. ((2T|) in the Virial expansion. The free 
energy per nucleon FjA is also shown for Hartree mean 
field calculations at intermediate densities, and from uni- 
form matter at high densities. The latter two have been 
obtained in our first paper Q- 

In most cases the transition (as density grows) is found 
at the density when the Hartree or uniform matter calcu- 
lation gives a lower free energy compared to Virial results. 
For matter at very low temperature (not greater than ~ 
1 MeV) and low proton fraction (not greater than ~ 0.1), 
some matching points are obtained at the density when 
the Virial gas and Hartree calculation give the closest free 
energy. The difference is below hundreds of KeV and is 
comparable to the mean deviation of binding energy for 
nuclei (- 600 KeV) in the FRDM mass table [Sj]. 



In each panel in Fig. [31 the vertical red dashed curves 
give the Virial gas - Hartree Wigner-Seitz cell transition 
densities, and the red solid curves give the transition den- 
sities to uniform matter. As temperature increases, the 
second transiton may happen at lower density than that 
of first transiton, which means the Hartree Wigner-Seitz 
cell is energetically unfavorable for all densities. This in- 
dicates the critical temperature for the nucleon liquid - 
Virial gas transition. We note here that the Virial gas 
may still contain a small amount of alpha particles and 
heavier nuclei even above this transition temperature. As 
seen from the figures, for matter at any temperature and 
proton fraction the Virial - mean field transition den- 
sities are larger than a few times 10"'' fm~^, which is 
about the neutrinosphere density. So in almost all re- 
gions around and below the neutrinosphere density, our 
EOS is represented by the Virial results, which include 
multiple nuclei. 

For matter at low density in the Virial gas phase, the 
free energy scales nicely with density. The reason is as 
follows. Low density matter is dominated by neutrons 
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and protons, and the interactions among them are not 
important due to the large particle spacing, so the kinetic 
energy dominates and scales as the logarithmic value of 
the density. As the temperature rises, the scaling behav- 
ior extends to higher densities until heavy nuclei appear. 



and electromagnetic and strong interaction corrections 
become important. Formation of nuclei greatly decreases 
the free energy, both in the Virial gas and Hartree mean 
field regimes. 
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FIG. 3: (color on line) Free energy per nucleon of nuclear matter at temperatures of T = 1 (a), 3.16 (b), 6.31 (c), and 10 (d) 
MeV. The proton fraction ranges from — 0.05 to 0.5. 



C. Mass fractions of species 

The Virial expansion gives the distribution of heavy 
nuclei (refer to Eq.(fT3l E])), where 8980 species of nu- 
clei are in thermal and chemical equilibrium with free 
neutrons, free protons and alpha particles. This is an 
improvement over the L-S EOS and S-S EOS that both 
used a single-nucleus representation. 

Figure m shows mass fractions of different nuclei (Z, A) 
for matter with T — 1 MeV, ub — 10"** fm"^, and Yp 
= 0.2. In upper panel different colors indicate the mass 



fraction using a Logj^g scale. In lower panel different lines 
indicate the contours of mass fraction in Logj^Q scale from 
-1 to -7. The total mass fraction of heavy nuclei is 56.5% 
(the rest is free neutrons). The majority of the mass frac- 
tions are centered around Z = 25 ^ 30, with a wide range 
of other nuclei with smaller abundances. In Figure [5] for 
a higher Yp=0.4, the total mass fraction of heavy nuclei 
is 99.1% (free neutrons and protons have 0.8% and 0.1% 
respectively). The majority mass fractions are centered 
around Z = 30 ^ 35. 

In Fig. [5] we show mass fractions of different nuclei for 
matter with T = 1 MeV, Yp = 0.4, and ub = 10"^ fm^^. 
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The total mass fraction of heavy nuclei is near one and 
the distribution is centered around Z = 35 and 50. The 
mass distribution of heavy nuclei in this high Yp case is a 
double-peaked Gaussian distribution, as shown in Fig. |S1 
where n{Z) is sum of the abundances of heavy nuclei with 
same proton number Z. At a higher T = 3.16 MeV shown 
in Fig. [71 the total mass fraction of heavy nuclei is 22.3% 
(free neutron, proton and alpha particles have 24.3%. 
6.3% and 47.1% respectively). Here the distribution of 
heavy nuclei is centered around a lower Z region and has 
multiple peaks in n{Z) as plotted in Fig. [S] It is also 
interesting to note that in this case there is odd-even 
effect in the value of Z for the abundances niZ). 

It is instructive to compare the composition of mat- 
ter in the Virial EOS with the existing EOS tables of 
Lattimer-Swesty and H Shen tt ai. The location of the 
neutrinosphere in a supernova is sensitive to the compo- 
sition of matter and is important for the emitted neu- 
trino spectra. Studies of collective flavor oscillations of 
neutrinos during their streaming outside of the neutri- 
nosphere have already indicated a sensitivity of neutrino 
oscillations to the emitted neutrino spectra from the neu- 
trinosphere (32| . Below we will compare some examples 
for the composition of matter around the neutrinosphere 
from the Virial EOS, the L-S EOS, and S-S EOS. 

Fig. [S] shows the mass fractions of neutrons, protons, 
alpha particles, and nuclei in matter at densities from 
10^^ to 10~^ fm~'^. The matter has a temperature of 
3.16 MeV and a proton fraction of 0.1 (a) or 0.3 (b), re- 
spectively. In panel (a) for a proton fraction of 0.1, free 
neutrons and protons dominate until the density reaches 
lO^'* fm^"^ in all three EoSs. Free nucleons dominate at 
low densities because of the high entropy. Above 10~^ 
fm^'^, alpha particles appear. The S-S EOS is close to 
our Virial results at densities below roughly 10"'^ fm~'^. 
The L-S EOS significantly underestimates Xq and this 
may be due to an error in the alpha particle binding en- 
ergy. Alpha particles have larger abundances and exist 
up to higher densities in our Virial EOS than in the other 
EOSs. This is partly because the attractive interactions 
between neutrons and alpha particles in the Virial expan- 
sion favors more alpha particles (17i | . Heavy nuclei begin 
to appear around 4x10"'' fm~^ in the L-S EOS, and 
at higher densities in the S-S EOS and our Virial EOS. 
Moreover, the L-S EOS predicts the largest abundance 
for heavy nuclei, while ours predicts the smallest abun- 
dance. Free neutrons have the largest abundance in our 
Virial EOS. This is due to the strong attractive interac- 
tion between neutrons in the Virial expansion which low- 
ers the energy and enhances the abundance of neutrons. 
Note in the Yp ~ 0.1 case the Virial-Hartree transition 
happens at 0.0158 fm^^^. The right panel of Fig. ^ for 
a different proton fraction of 0.3, has similar characteris- 
tics. However here, alpha particles and heavy nuclei have 
much larger abundances than for the Yp = 0.1 case, since 
a higher proton fraction favors formation of nuclei. In 
this Yp — 0.3 case, the transition density from Virial gas 
to Hartree mean field calculations occurs at 6.3 x 10"'^ 



fm ^ as indicated by the dotted line in the figure. 

In future work 33] we will carefully interpolate the free 
energy results of this paper in a thermodynamically con- 
sistent way in order to accurately determine derivatives 
of the free energy with respect to temperature or den- 
sity. From these derivatives we will calculate a number 
of additional quantities such as the pressure or entropy. 



V. SUMMARY AND OUTLOOK 

In this paper we present our model for nuclear matter 
at subnuclear density, the Virial expansion for a nonideal 
gas consisting of neutrons, protons, alpha particles and 
thousands of heavy nuclei. We include second order virial 
corrections for light elements A < 4, nuclear partition 
functions for heavy nuclei, and Coulomb corrections. At 
very low density, the Virial expansion reduces to nuclear 
statistical equilibrium. We calculate the free energy, and 
tabulate the resulting EOS at over 73,000 grid points 
in the temperature range T = 0.158 to 15.8 MeV, the 
density range ub = 10~^ to 0.1 fm"^, and the proton 
fraction range Yp = 0.05 to 0.56. These calculations 
took over 1000 CPU days. The above parameter space 
is complementary to that of the relativistic mean field 
results in our previous paper 0]. 

The treatment of Coulomb corrections in Wigner-Seitz 
approximation agrees reasonably with an analytical clus- 
ter expansion. Our results do not appear to be very sen- 
sitive to the mass table employed, or to the form of the 
partition function for heavy nuclei. However, results at 
higher densities for mean field calculations are sensitive 
to the interaction employed. In the future, we intend to 
match the Virial results of this paper to mean field cal- 
culations using a number of different interactions, see for 
example [3^ . 

Our EOS includes broad distributions of nuclei, that 
are not included in two commonly used EOS tables. 
These distributions of nuclei make the composition of nu- 
clear matter in our EOS different from L-S and S-S EOS 
tables. These differences may be important for neutrino 
interactions. 

This paper provides the second part of our results for 
a complete EOS table that will cover a broad range of 
temperatures, densities, and proton fractions for use in 
supernova and neutron star merger simulations. In the 
future we will use a thermodynamically consistent inter- 
polation scheme to match the Virial EOS in this paper 
and the RMF EOS from our previous paper, or similar 
RMF models, and generate a complete EOS table (ssj . 
Our Virial EOS is exact in the low density limit. Finally 
the Virial expansion also provides a suitable framework 
for future work that includes other modern mass tables, 
such as HFB14 [31] or KTUY05 |35|]. 
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FIG. 4: (color on line) Mass fraction of nuclei in the nuclear chart for matter at T = 1 MeV, ub ~ 10~^ fm^'', and Yp — 0.2. 
Upper panel: different colors indicate the mass fraction using a Logj^Q scale. Lower panel: different lines indicate the contours 
of mass fraction in Logj^Q scale from -1 to -7. 
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FIG. 5: (color on line) Mass fraction of nuclei in the nuclear chart for matter at T = 1 MeV, ns — 10~^ fm^'', and Yp — 0.4. 
Upper panel: different colors indicate the mass fraction using a Logj^Q scale. Lower panel: different lines indicate the contours 
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FIG. 7: (color on line) Mass fraction of nuclei in the nuclear chart for matter at T = 3.16 MeV, ub ~ 10~^ fm~'^, and Yp 
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FIG. 8: (color on line) Mass fractions of nuclei n{Z) for nuclei proton number Z, for the cases in Figs.|4l[5l[6]and[71 The triplet 
of values in the parenthesis are (r/[MeV], Yp, nB/[im-^]). 




FIG. 9: (color on line) Mass fractions of matter at T = 3.16 MeV, Yp = 0.1 (a) and 0.3 (b). The solid curves show our Virial 
EOS results, while the dashed lines show Lattimer-Swesty EOS results and the H Shen et al. EOS results are shown by the 
dot dashed curves. For Yp = 0.3, the dotted line indicates the transition density between Virial EOS and Hartree mean field 
results. For Yp — 0.1, the transition density between Virial EOS and Hartree is 0.0158 fm"''. 



